Kosterlitz-Thouless transition in disordered two-dimensional topological insulators 
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The disorder-driven metal-insulator transition in the quantum spin Hall systems is studied by 
scaling analysis of the Thouless conductance g. Below a critical disorder strength, the conductance 
is independent of the sample size M, an indication of critically delocalized electron states. The 
calculated beta function ft = din g/ din M indicates that the metal-insulator transition is Kosterlitz- 
Thouless (KT) type, which is characterized by bounding and unbounding of vortex-antivortex pairs 
of the local currents. The KT like metal-insulator transition is a basic characteristic of the quantum 
spin Hall state, being independent of the time-reversal symmetry. 
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The Kosterlitz-Thouless (KT) transition, which was 
first proposed by Kosterlitz and Thouless in two- 
dimensional (2D) XY model, [1] is an important phe- 
nomenon in 2D systems, such as 2D magnets, supercon- 
ductors and superfluids in thin films. It is a typical topo- 
logical phase transition, which has been understood as 
bounding and unbounding of vortex-antivortex pairs. [I], 
2] In recent years, the KT like transition was also found 
in 2D disordered electronic models, including the model 
with random magnetic fluxes, [3, [ij and the graphene 
model with long-range impurity scattering potential. [f| 
These systems show localization-delocalization transi- 
tions with decreasing disorder strength. On the metallic 
side, the Thouless conductance g is independent of the 
system size M, leading to a vanishing beta function [1,0] 
defined as f3 = dlng/dlnM, which is an essential charac- 
teristic of the KT transition. From the classification 
scheme based on the universality classes, @ the former 
model [1, 0] belongs to the unitary class, and the exis- 
tence of delocalized states is due to the random phases in 
the hopping integrals, which result in the Chern number 
fluctuations. [|[ The latter model @ belongs to the or- 
thogonal class, and the electron derealization arises from 
the nonzero Berry phases at two Dirac points, which can- 
not annihilate each other when the impurity scattering 
is long-range correlated. 

A topological insulator [l(| [H| is a band insulator with 
gapless edge states or surface states traversing the band 
gap, which originate from the nontrivial band topology. 
The gapless nature of the edge states or surface states is 
protected by the time-reversal symmetry. The 2D ver- 
sion of the topological insulators, which is also called the 
quantum spin Hall (QSH) state, was predicted by Kane 
and Mele [12j | to exist in a graphene model with intrinsic 
spin-orbit coupling. Later it was experimentally realized 
in mercury telluride (HgTe) quantum wells. [l3[ The ro- 
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bustness of the QSH state is associated with new types of 
topological invariants of the bulk energy bands. [l4Hl7| 

Since topological invariants are carried by bulk ex- 
tended states, there have being strong interest in the 
derealization problem in the QSH systems. Onoda et 
al. [H[ studied the Kane- Mele model using the trans- 
fer matrix method and showed the existence of extended 
states in a finite energy region in the presence of time re- 
versal symmetry. Based upon calculations of the critical 
exponents, they suggested that the QSH state belongs 
to a new universality class different from the ordinary 
symplectic universality class. [l!| Consistent result was 
obtained by Prodan using the level statistics analysis. [2(| 
Employing a continuum model, Mong et al. (2l| showed 
that the surface states of a weak topological insulator 
with even number of Dirac cores are delocalized in the 
absence of a mass term, and localization occurs when the 
mass term and strong disorder are present. This conclu- 
sion may also be applicable to the QSH systems. Xu et 
al. [13] demonstrated numerically that the metallic state 
in the QSH systems remains to be robust even when the 
time-reversal symmetry is broken, in glaring contrast to 
trivial unitary class, where all electron states are local- 
ized. This result was attributed to the fact that the non- 
trivial band topology, as described by the spin Chern 
numbers, is intact when the time-reversal symmetry is 
broken. 17, 23| The existence of delocalized states in the 
QSH systems is affirmative, but the basic characteristics 
of the metal-insulator transition have not been well es- 
tablished. 

In this Letter, from calculations of the Thouless con- 
ductance in the Kane-Mele model, it is shown that the 
disorder-driven metal-insulator transition in the QSH 
system is a KT phase transition. Below a critical disor- 
der strength, the Thouless conductance g is found to be 
independent of the sample size, an indication of critically 
delocalized electron states. The beta function obtained 
from the universal scaling function of the conductance, 
vanishes in the metallic phase. Through mapping the lo- 
cal currents onto the local spins in the 2D XY model, 
we further show that the KT transition is characterized 
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by bounding and unbounding of vortex-antivortex pairs. 
The KT like metal-insulator transition is a fundamental 
characteristic of the QSH state, allowing to be exam- 
ined experimentally, which neither depends on the time- 
reversal symmetry, nor requires any long-range correla- 
tions of the disorder. 

We start from the Kane-Mele model defined on a 2D 
honeycomb lattice, [IH, [ll| [24j with the Hamiltonian 
given by 

H = -t^2c\cj + ^teV so 4° ■ (dkj x dik)cj 

(ij) V ««» 



iV a ^ ' (<* x dij)cj + ^2 WiCi'Ci . (1) 



(ij) 



Here, the first term is the usual nearest neighbor hopping 
term with q = (cjj., c]^) as the electron creation operator 
on site i. The second term is the intrinsic spin-orbit cou- 
pling with coupling strength V so , where a are the Pauli 
matrices, i and j arc two next nearest neighbor sites, k 
is their unique common nearest neighbor, and vector dik 
points from k to i. The third term stands for the Rashba 
spin-orbit coupling with coupling strength V R . Wi is a 
random on-site energy potential uniformly distributed 
between [—W/2,W/2], which accounts for nonmagnetic 
disorder. For convenience, we will set h, t and the dis- 
tance between the nearest neighbor sites all to be unity. 
In the following calculations, V^ = V R = 0.1 is fixed, so 
that the model is inside the QSH phase at W = 0. 

The localization property of an infinite 2D system can 
be extracted from the scaling analysis of the dimension- 
less Thouless conductance calculated for finite-size sam- 
ples. [1, Q We consider a rectangular sample with two 
zigzag sides and two armchair sides. Each zigzag side 
has M atoms and each armchair side has 2M atoms, so 
that the sample has totally M x M atoms. The rectangu- 
lar sample is connected to two semi-infinite leads of the 
same width with armchair interfaces. At zero tempera- 
ture, the two-terminal dimensionless Landauer conduc- 
tance Or is calculated from the Landauer-Biittiker for- 
mula [25] g L = Tr[r L G"T R G a ], where G r ^{E) is the 
retarded (advanced) Green's function and T UR) {E) is 
the coupling matrix between the sample and the left 
(right) lead. The Thouless conductance g is related 
to the Landauer conductance g^ through the relation 
1/.9 = VSl - I /N a , @, 1H where 1/N a is the contact 
resistance between the leads and the sample with N c as 
the number of propagating channels in the leads at the 
Fermi energy E. 

Figure 1 shows the size dependence of the calculated 
conductance g in logarithmic scale for different disorder 
strengths W at E = —1 and —2. Periodic boundary 
conditions are employed in the transverse direction, so 
that no contribution from edge states is involved. At 
given values of W and M, the conductance is averaged 
over 200 to 500 random disorder configurations. Sym- 
bols with different shapes are used to distinguish data 
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FIG. 1: (Color online) Conductance as a function of system 
size M at E = — 1 (a) and E — — 2 (b). Symbols with different 
colors and shapes correspond to different disorder strengths 
W. 



for different W . It is found that the conductance behaves 
differently with changing sample size in weak and strong 
disorder regions. For weak disorder, g neither decreases 
nor increases with increasing sample size. As a result, 
the conductance is expected to remain constant in the 
thermodynamic limit, suggesting that the electron states 
at the Fermi energy are critically delocalized. For strong 
disorder, g decreases with increasing sample size, corre- 
sponding to localized states. A metal-insulator transition 
occurs at an intermediate critical disorder strength. 

To further investigate the nature of the metal-insulator 
phase transition, we perform scaling analysis on the 
data shown in Fig. 1 and calculate the beta function 
/3 = d\ng/d\nM. According to the one-parameter scal- 
ing theory, [10] 

the localization length in the thermo- 
dynamic limit can be extracted by fitting all the data 
on the insulator side in Fig. 1 into a universal function 
g = f(M/£), where the fitting parameter £ yields the 
localization length in the thermodynamic limit. The ob- 
tained localization lengths for E = — 1 and —2 are plotted 
in Figs. 2a and 2b, respectively, and the universal scal- 
ing function g = /(M/£) is shown in Fig. 2c. In the 
above scaling procedure, only the relative magnitudes of 
the localization lengths £ can be determined. Therefore, 
we set £ at W = 6.0 to be unity, and the localization 
lengths for other W are the relative values with respect 
to it. It is found that the localization length £ as a func- 
tion of disorder strength W shown in Figs. 2(a,b) can 
be fitted with an exponential function £ oc e a /VW-w c 
with a and W c as two fitting parameters. This behavior 
indicates an exponential decay of the localization length 
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FIG. 2: (Color online) (a,b) Calculated localization length £ 
(symbols), fitted with an exponential function (solid lines), 
(c) Universal scaling function g — /(M/f), where solid and 
open symbols are conductance at E = — 1 and —2, respec- 
tively, and the solid line represents a smooth functional fit to 
the data, (d) Calculated beta function. 



with increasing the disorder strength in the insulating 
phase, which is typical of a KT phase transition, [l|, [2| 
The critical disorder strength for such a metal-insulator 
transition is obtained as W c — 3.12 at E = — 1 and 2.54 
at E = —2, respectively. To facilitate the calculation of 
the beta function /3 = <iln g/d\n M, we fit the data in 
Fig. 2c with a smooth function, as indicated by the solid 
line, which is explained below. In the region of small 
ln(M/£), the data is fitted with g = g c e~ ( - M ^'> 1/ \ which 
is consistent with the general behavior of the beta func- 
tion, (3 — i ln(gr/<7 c ), in the quantum critical region. (27j 
On the other hand, a polynomial fit is used in the region 
of relatively large ln(M/£). At the boundary between 
the two regions, In g and its derivative f3 are required to 
be continuous. 

After obtaining the smooth fitting function describ- 
ing the scaling relation between In g and InM, we 
take its derivative to calculate the beta function /? = 
d\ng/dhiM, as plotted in Fig. 2d. Here, we note that 
for weak disorder {W < W c ), where the conductance 
is nearly independent of the sample size, the localization 
length diverges and cannot be treated with the same scal- 
ing procedure as above. Instead, we fit the logarithmic 
conductance In g as a function of In M with straight lines, 
and take their slopes as the values of j3, which are shown 
by the diamonds in Fig. 2d. We see from Fig. 2d that (3 
increases with lng, and vanishes on average above a crit- 
ical point (g > g c ~ 1.3) with negligible fluctuations, in- 
dicating clearly that the disorder-driven metal-insulator 
transition is a KT transition, as observed in other disor- 
dered 2D electron systems. The behavior of the beta 
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FIG. 3: (Color online) Typical distribution of vortices (red 
positive numbers) and antivortices (blue negative numbers) 
on the two sides of the KT like transition in a sample with 
40 x 40 sites at E — — 2 for (a) W = 1.5t (metallic phase) 
and (b) W = 3.5t (insulator phase). For clarity, only a part 
of the sample is shown. 



function also suggests that a minimum (critical) conduc- 
tivity, a c = g c {e 2 /h)/V3 ~ 0.75(e 2 //i) (l/%/3 being the 
aspect ratio), exists for the bulk QSH system. 

An important characteristic of the traditional KT tran- 
sition in the XY model is the bounding of vortex- 
antivortex pairs in the ordered phase near the transition 
point and unbounding in the disordered phase, [l], Q For 
2D electron systems, Zhang et al. [5[ proposed to map 
the local currents onto the local spins in the XY model. 
The bond current vector between sites m and n can be 
calculated by @, H, EH 



l, m.(T2— *na\ 



-2e f dE 



l^wi.muj^)]) ( 2 ) 



where o\ and oi represent the spin index and 



, (E) is the matrix element of the lesser Green's 



function. The lesser Green's function is given by 

G<(E) =G r (E)[iY / r a (E)f a (E)]G a (E) , (3) 



where f a (E) = fo(E + eV a ) is the Fermi distribution 
function in lead a with V a as the electrical potential in 
the lead. The local current vector defined on site n is 
in = Y, m <7 2 <Ji Wa-Mwn @ where the vectorial summa- 
tion is taken over all the nearest and next nearest neigh- 
bors m of site n. Along a closed path, the polar angle 
9 n of i n is considered to change continuously with chang- 
ing coordinates, and a topological charge is defined as 
c= <f V9 ■ dl. [j| A nonzero number of c indicates the 
appearance of a vertex (c > 0) or an antivortex (c < 0). 
Figure 3 shows typical distributions of the topological 
charges on the two sides of the metal-insulator transi- 
tion. The arrows in the figure indicate only the direction 
of the local currents, without showing the magnitude for 
the purpose of visualization. In the delocalized phase 
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FIG. 4: (Color online) (a-d) Conductance as a function of 
sample size with Fermi energy E, in the presence of a vertical 
Zeeman field hz or a horizontal Zeeman field h x . (e) Universal 
scaling function g — /(M/£), where the solid line represents a 
smooth functional fit to the data, (f) The corresponding beta 
functions for the time-reversal symmetry broken QSH system 
(solid line) and time-reversal symmetry invariant QSH system 
(dotted line). 

(Fig. 3a), a few vortices (c > 0) and antivortices (c < 0) 
are excited in pairs and bounded together, corresponding 
to the ordered phase at low temperature in the 2D XY 
model. In the localized phase (Fig. 3b), a large number 
of vortices and antivortices appear and are mostly un- 
bounded, corresponding to the disordered phase at high 
temperature. 

Finally, we study the effect of breaking the time- 
reversal symmetry on the KT phase transition, by in- 
cluding a new term h ■ a into Hamiltonian ([1]), where h 
stands for a uniform Zeeman field. The calculated con- 
ductances at E = — 1 and —2 as functions of sample 
size M for different disorder strengths are shown in Fig. 
4(a,b) for a vertical Zeeman field h = (0, 0, h z = 0.1) and 
in Fig. 4(c,d) for a horizontal field h = (h x — 0.1,0,0). 
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The parameters are chosen in such a manner that the sys- 
tem is in the time-reversal symmetry broken QSH phase 
in the clean limit. [22|, [23[ The conductance shown in Fig. 
4 (a-d) displays similar size dependence to that in Fig. 1. 
For weak disorder, the conductance is nearly independent 
of sample size, indicating the existence of the critically 
delocalized states. For strong disorder, the conductance 
decreases with sample size, indicating that the electron 
states are localized. It is found that all the data on the 
insulating side in Fig. 4(a-d) can also be fitted with a 
universal scaling function g — /(M/£), as shown in Fig. 
4e. In Fig. 4f, the beta function (solid line) for the time- 
reversal symmetry broken QSH system is obtained in the 
same way as in Fig. 3. As expected, the beta function 
vanishes above a critical conductance g c , exhibiting the 
essential feature of the KT like transition. For compar- 
ison, /3 for the time-reversal invariant QSH system ob- 
tained in Fig. 2d is shown by the dotted line in Fig. 4f. 
Both curves exhibit the same characteristic of the KT 
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In summary, the metal-insulator transition in the QSH 
systems is studied based upon the scaling analysis of 
the Thouless conductance. The disorder-driven metal- 
insulator transition is found to be a KT like transition, 
characterized by bounding and unbounding of vortex- 
antivortex pairs of local electrical currents. The beta 
functions for the QSH systems with and without the 
time-reversal symmetry are close to each other on the 
metallic side, but somewhat different on the insulating 
side, the latter being attributed to the fact that they be- 
long to different symmetry classes. Our work established 
a fundamental property of the QSH state, which can be 
observed experimentally. 
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